%% Change_thickness_CF
clear all; clc;
E1=350*10^9;  E2=7*10^9; E=62*10^9;  E3=300;
v12=0.33; v21=v12*E2/E1; v=0.34;
d31=-320*10^(-12);  d32=d31;  d33=0;
L=0.01;   W=0.002;    H=127/10^6;
rho=7800;
m=L*W*H*rho;
%% 
gamma=45;   % designated ply angle
%% 
gamma_1=-gamma*pi/180;
m1=cos(gamma_1);
n1=sin(gamma_1);
gamma_2=gamma*pi/180;
m2=cos(gamma_2);
n2=sin(gamma_2);
T1=[m1^2 n1^2 2*m1*n1;n1^2 m1^2 -2*m1*n1;-m1*n1 m1*n1 m1^2-n1^2];
T2=[m2^2 n2^2 2*m2*n2;n2^2 m2^2 -2*m2*n2;-m2*n2 m2*n2 m2^2-n2^2];
O11=E/(1-v^2);
O12=E*v/(1-v^2);
O22=O11;
O66=E/(2+2*v);
Q11=E1/(1-v12*v21);
Q12=v12*E2/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=350*10^9;
Q_bar=[Q11 Q12 0;Q12 Q22 0;0 0 Q66];
q1=inv(T1)*Q_bar*(inv(T1))';
q2=inv(T2)*Q_bar*(inv(T2))';
%% Change the thickness of passive layer of Carbon fiber
for t=1:200            % t=1*10e-6:200*10e-6
z1=-(127/10^6)/2;
z2=-z1;
z3=z2+t/10^6;
z0=-z3;
A11=q1(1,1)*(z1-z0)+q2(1,1)*(z3-z2)+O11*(z2-z1);
A12=q1(1,2)*(z1-z0)+q2(1,2)*(z3-z2)+O12*(z2-z1);
A22=q1(2,2)*(z1-z0)+q2(2,2)*(z3-z2)+O22*(z2-z1);
A21=A12;
A31=q1(3,1)*(z1-z0)+q2(3,1)*(z3-z2);
A13=A31;
A32=q1(3,2)*(z1-z0)+q2(3,2)*(z3-z2);
A23=A32;
A33=q1(3,3)*(z1-z0)+q2(3,3)*(z3-z2)+O66*(z2-z1);
B11=0.5*q1(1,1)*(z1^2-z0^2)+0.5*q2(1,1)*(z3^2-z2^2)+0.5*O11*(z2^2-z1^2);
B12=0.5*q1(1,2)*(z1^2-z0^2)+0.5*q2(1,2)*(z3^2-z2^2)+0.5*O12*(z2^2-z1^2);
B22=0.5*q1(2,2)*(z1^2-z0^2)+0.5*q2(2,2)*(z3^2-z2^2)+0.5*O22*(z2^2-z1^2);
B23=0.5*q1(2,3)*(z1^2-z0^2)+0.5*q2(2,3)*(z3^2-z2^2);
B13=0.5*q1(1,3)*(z1^2-z0^2)+0.5*q2(1,3)*(z3^2-z2^2);
B33=0.5*q1(3,3)*(z1^2-z0^2)+0.5*q2(3,3)*(z3^2-z2^2)+0.5*O66*(z2^2-z1^2);
B21=B12;
B31=B13;
B32=B23;
D11=1/3*q1(1,1)*(z1^3-z0^3)+1/3*q2(1,1)*(z3^3-z2^3)+1/3*O11*(z2^3-z1^3);
D12=1/3*q1(1,2)*(z1^3-z0^3)+1/3*q2(1,2)*(z3^3-z2^3)+1/3*O12*(z2^3-z1^3);
D22=1/3*q1(2,2)*(z1^3-z0^3)+1/3*q2(2,2)*(z3^3-z2^3)+1/3*O22*(z2^3-z1^3);
D33=1/3*q1(3,3)*(z1^3-z0^3)+1/3*q2(3,3)*(z3^3-z2^3)+1/3*O66*(z2^3-z1^3);
D13=1/3*q1(1,3)*(z1^3-z0^3)+1/3*q2(1,3)*(z3^3-z2^3);
D23=1/3*q1(2,3)*(z1^3-z0^3)+1/3*q2(2,3)*(z3^3-z2^3);
D21=D12;
D31=D13;
D32=D23;
A=[A11 A12 A13;A21 A22 A23;A31 A32 A33];
B=[B11 B12 B13;B21 B22 B23;B31 B32 B33];
D=[D11 D12 D13;D21 D22 D23;D31 D32 D33];
C=inv([A B;B D]);
N1=(O11*d31+O12*d32)*E3*(z2-z1);
N2=(O12*d31+O22*d32)*E3*(z2-z1);
N3=(O66*d33)*E3*(z2-z1);

kxy=C(6,1)*N1+C(6,2)*N2+C(6,3)*N3;
theta(t)=atan(kxy*L)*180/pi;
tau(t)=W*(C(6,1)*N1+C(6,2)*N2+C(6,3)*N3)/C(6,6);
Ed(t)=0.5*tau(t)*theta(t)/m;
end
t=1:200;
figure(1)
plot(t,theta)
figure(2)
plot(t, tau)
figure(3)
plot(t, Ed)


